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ABSTRACT 



Acoustic ray paths in the ocean are known to exhibit sig- 
nificant horizontal deflections after repeated reflections 
from the bottom. The effect may be quantitatively and quali- 
tatively observed through a ray trace model which permits a 
change in direction of the vertical plane of propagation as a 
function of bottom slope and grazing angle. ECTRACE is a 
family of computer programs which traces a bundle of rays in 
three dimensions and utilizes bottom depth values as a portion 
of its input. Included are related programs which develop sea 
bed models from digitized bathymetry data or synthetic bathymetry 
functions . 
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I . INTRODUCTION 



ECTRACE is a FORTRAN computer program designed for the 
investigation of bathymetric effects on horizontal acoustic 
ray tracing. Its supporting computer programs include routines 
for constructing a grid model of the sea bed from contour data, 
ray tracing routines, and associated printed and plotted output 
The programs were written for IBM 360/370 series systems with 
a Versatec plotter and related software. The significant fea- 
ture of the tracing algorithm is its ability to shift the direc 
tion of the vertical plane of propagation after a bottom 
reflection. A sea bed modeled by triangular facets fitted 
through three points of a grid cell makes this feature possible 
with a minimum effect on computation time. 

Several models for tracing acoustic rays in the ocean are 
used in the Navy today. These models are based on a simplified 
environment because of practical limitations on data and com- 
puter run time, and they adequately trace a bundle of rays in 
a single vertical plane. Some enable intensity calculations 
at a given distance from the source. TRIMAIN [1], for example, 
is a model which features range as well as depth dependence 
of sound speed and accounts for some interaction with irregu- 
lar bathymetry. These models have a large number of applica- 
tions in deep ocean cases involving ducts, channels, and 
convergence zones, and bottom reflections where the plane of 
the ray path does not exhibit horizontal deflection. 
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Actual cases of long range propagation over irregular bot- 
toms show any single ray path rarely remains within a single 
vertical plane. Any application of ray theory to propagation 
involving bottom interaction must account for a horizontal 
deflection caused by ray heading changes after repeated reflec- 
tion from a tilted bottom plane. This effect is most pro- 
nounced for long range paths over smoothly sloping surfaces, 
but any valid ray approximation can be seen to undergo signif- 
cant heading changes after only a few reflections from an 
irregular or undulating bottom. 

The behavior of sound propagation over troughs, ridges, 
wedges, and seamounts has been studied analytically by Har- 
rison [2] using normal mode theory and ray invariants. Since 
actual ocean bottom topography over large area defies prac- 
tical analytical description, except in the stochastic sense, 
DeWitt [3] developed numerical techniques for three-dimensional 
ray tracing in which the bottom consists of triangular facets 
generated numerically from actual or interpolated bathymetry 
data. This approximation scheme greatly simplified the prob- 
lem of calculating the three-dimensional ray path parameters 
after a bottom reflection and provides a means of constructing 
reasonably accurate ray traces (for a verification, see Appen- 
dix F) . 

The ECTRACE modeling program is an adaptation of the 
Dewitt technique for use on an IBM 360/370 series computer. 

The initial programs may be used to generate synthetic bot- 
tom topographical grids or create approximations to real 
ocean bottom grids from a bathymetric data source file. The 
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topographical features of the generated grid can be depicted 
by a contour plot and a perspective surface plot for compari- 
son with the original contours . The primary program is the 
ray tracing program which accepts the generated bottom grid 
and performs stepwise ray position calculations from user- 
selected initial heading and elevation angles. The output 
is the printed history of each ray, and a two-dimensional plot 
of the ray paths projected on a vertical plane intersecting 
the bottom profile along the axis of the ray fan. Addition- 
ally, a horizontal plane projection is plotted which may be 
overlaid on the appropriate contour plot to correlate the 
horizontal curvature effect with grid bathymetry. Punched 
card output from several ray tracing runs may be combined in 
one independent routine for a comprehensive visual inspection 
of acoustic shadow zones in the horizontal plane. 

It is very important to note that ECTRACE permits refrac- 
tion caused by sound velocity gradients in the vertical plane 
only and that the sound velocity in the horizontal plane is 
assumed constant within a well defined water mass. All hori- 
zontal bending effects revealed by ECTRACE are the result of 
bottom reflection only. 

With the ability to convert actual bathymetry data to a 
matrix of depth values, discrete faceted approximations of a 
selected oceanic region can be used for generating ray traces 
of operational significance. In addition, qualitative studies 
of horizontal ray deflections over idealized bottom configura- 
tions are enhanced by the computer-generated plots of ray paths 
over easily modifiable synthetic bathymetry. ECTRACE and the 
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augmenting programs are envisioned to serve as a basis for 
future development of three-dimensional ray modeling and inves- 
tigation efforts. 
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II. DESCRIPTION OF ECTRACE AND PERIPHERAL PROGRAMS 



This chapter describes each of the programs, their pri- 
mary subroutines, and their interrelationship in producing a 
ray trace which may be used to investigate the bathymetric 
effects on horizontal ray curvature. They are presented in 
their intended order of use. Flow Diagram (Fig. 1) and Table 
I are provided to assist in understanding program relationships. 

A. GENBOT 

This program converts contour data into a matrix of depth 
values for storage on a permanent device. It also makes a 
two dimensional contour plot of the input data and super- 
imposes reference latitude and longitude lines (Fig. 2). 

The contour plot can be placed underneath the horizontal ray 
trace plot of ECTRACE or ECCOM on a light table to aid in 
the visual investigation of the ray curvature (see Figs. 3 
and 4, and Figs. 5, 6, and 7). 

The bathymetry contour data used in this research exists 
as a sequence of data sets representing ten by ten degree 
regions of the North Atlantic Ocean and were obtained exclu- 
sively from a tape compiled in a joint project of the Naval 
Research Laboratory and the Mobil Oil Company. In their 
original form the data were sufficient only for producing a 
computer plot of region contours, thereby reproducing por- 
tions of a chart titled "Bathymetry of the Norwegian-Green- 
land and Western Barents Sea" [4] (see Fig. 8). Their use 
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in GENBOT required us to insert the depth values of the con- 
tour lines and transfer the results to a separate tape (de- 
scribed in Appendix A) . 

The following subroutines are called by GENBOT: 

1) CORNER calculates boundary latitudes and longitudes for 
screening only those data needed for interpolation. 

2) GEOPLT superimposes reference latitude and longitude 
lines on the subregion contours. 

3) GEODST calculates the forward or inverse solutions of 
the geodetic triangle to accomplish the above. 

4) RAIN 1, RAIN 2, and RAIN 3 [5] were obtained from the 
SSP3 program library at NPS and are used to perform 
interpolations on data points for construction of the 
depth matrix. 

B. SYNGEN 

This program produces a synthetic bathymetry grid for 
storage on a permanent device. It may be used to generate 
the following types of ideal sea bed configurations: wedge, 

trough, ridge, conical seamount (Fig. 9), sinusoidal undula- 
tion, or parabolic basin (Fig. 10) . These grids can be used 
in the ECTRACE program as a depth matrix in the same manner 
as the grids of actual bathymetry produced by GENBOT. SYNGEN 
is discussed in detail in Appendix B. 

C. G3DP 

Program G3PD produces a perspective of any portion of the 
generated depth matrix through use of the NPS system subroutine 
CONTUR [4], Figures 9, 10, and 11 are examples. G3DP uses 
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job control commands to link with the GRDSCT subroutine of the 
ECTRACE load module to extract the desired area (hereafter 
called the working matrix) to be plotted. Appendix C describes 
the application of G3DP in detail. 

D. GRDCHK 

This program is used to check the validity of the depth 
matrix generated by the programs GENBOT and SYNGEN. It produces 
a contour plot and a vertical plot of user-selected rows of the 
generated matrix. By comparing consecutive rows the existence 
and location of anomalous data points may be determined (Figs. 

4 and 12) . It contains the previously mentioned subroutine 
GEODST, CORNER, and a version of GEOPLT (called GEOPLR) to draw 
reference geographic lines for comparison of the two contour 
plots (GENBOT vs GRDCHK) . Refer to Appendix D for additional 
discussion of GRDCHK. 

E . ECTRACE 

The ECTRACE program is a package of subroutines which 
exist in load module form on a permanent storage device in 
the NPS computer system user library. Its use requires a 
user supplied calling program which calls subroutine TRACER 
and contains the necessary job control language (JCL) to link 
with the module. These procedures are explained in more de- 
tail in Chapter IV of this report. 

ECTRACE traces rays in three dimensions. The sound speed 
field is assumed to vary piecewise-linearly with depth only, 
yet provisions are made to permit simulation of up the three 
distinct water masses separated discontinuously by vertical 
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fronts. The bathymetry structure, generated by one of the 
two grid-making programs, is a three dimensional surface 
modeled by discrete triangular facets fitted through cells 
of adjacent depth values. Rays are traced as piecewise arc 
segments each with a radius of curvature dependent on the 
vertically-structured sound speed gradient. 

When the rays undergo a bottom bounce, specular reflection 
is assumed. Bottom loss values can be calculated from the 
subprogram FBLOSS (described later) , a standard Navy model 
such as the FACT bottom-loss routine, or a user supplied 
routine. The trace history supplied by ECTRACE includes the 
location and characteristics of each ray reversal (surface 
reflections, bottom bounces, and refractive turning points). 
The printed data lists the new ray parameters determined 
after the reversal had occurred. Chapter V explains the 
printed output in more detail. 

The ECTRACE plot product includes a plane view of the 
ray paths projected on a vertical plane, a plot of sound 
speed profiles, and a horizontal (plane) view of the ray 
paths (Figs. 5 and 13). The horizontal view can be overlaid 
on the appropriate contour plot generated by GENBOT or GRDCHK, 
as in Figs. 12 and 14. 

The following subroutines are used by ECTRACE: 

1) Tracer is the primary tracing subroutine which calls 
all other subroutines of the ECTRACE program and is 
capable of processing a bundle of rays separated into 
ray fans. A ray fan is defined here as a number of 
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rays (user determined) each with a different initial 
heading measured clockwise from grid north, but all 
with the same initial elevation angle. A positive 
elevation angle describes a ray which is pointed toward 
increasing depth. 

2) GRDSCT reads into core a fixed-dimension section of a 
depth matrix file of arbitrary size. The ray trace is 
restricted to the boundaries of the core-loaded working 
matrix. A ray which departs the working matrix boundaries 
is terminated and program control is passed to the next 
ray's trace. 

3) IDTSUB identifies the three matrix subscript pairs whose 
depth values define a plane triangular bottom facet 
underneath the ray segment head. 

4) DEEP calculates the depth at any horizontal position by 
solving for the z-coordinate on a plane fitted through 
the three depth matrix values. 

5) CONTAC calculates the three-dimensional coordinates of 
the intersection of the ray path and triangular bottom 
facet . 

6) BOANG calculates the grazing angle, elevation angle, and 
new heading of a bottom reflected ray. 

7) FBLOSS returns bottom loss values in dB as a function 
of grazing angle for a bottom reflected ray. These 
values are based on NRL geophysical survey data of the 
Greenland-Norwegian Sea and Iceland environment collected 
in the mid-1970's and have been widely averaged for 
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convenience. The user has the option of substituting a 
more specialized function if bottom loss values are 
critical . 

8) IDPROF determines the water mass, and hence the sound 
speed profile, affecting the ray with each segment 
iteration . 

9) NUPROF initializes calculations for new water mass param- 
eters once the ray has passed the boundary. 

10) CHNLIM determines in advance the refractive sound channel 
limits for each ray from its invariant. This is used to 
test a ray at its turn-around point and to reveal its 
type (surface reflected, refracted-surface reflected, 
bottom ref lected-surface reflected) in the printed 
trace history. 

11) ANGPRT is used to print a summary of selected bottom 
contact in parameters for each ray fan. 

12) The following subroutines each have only one calling 
statement in TRACER. The subroutines they in turn invoke 
are dependent upon the computer graphic system installed 
and would require internal modification if exported to 

a system without Versatec software. For additional 
information NPS users should consult Ref. 7. 

a) BGNPLT initializes the plotter and draws all borders, 
titles, and axes. 

b) BDTPLT draws a vertical profile of the sea bed cen- 
tered along the mean heading of the ray bundle, as 
shown in the upper plot of Figure 13. 
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c) RNGPLT traces all ray segment projections on the 
vertical plot. 

d) T2DPLT plots the horizontal track of each ray. Each 
ray fan uses one of thirteen symbols to indicate the 
bottom bounce positions of each ray. Figure 5 shows 
some of the available symbols. 

e) SSPPLT draws the sound speed profiles. 

f) ENDPLT draws the horizontal plot legend and termi- 
nates all plotting. 

F . ECCOM 

Program ECCOM uses the optionally produced punched card 
output of one or more ECTRACE jobs to draw a composite hori- 
zontal plot of several ray bundles, enabling a comprehensive 
pictorial study of shadow zones and acoustic convergence in 
the horizontal plane. The composite horizontal plot is de- 
signed to portray the rays which emanate from the same source 
or from multiple sources along a desired track. The ray num- 
bers assigned and punched for each ray on each ECTRACE run 
assist the user in singling out rays of interest (or non- 
interest) before combination. Figures 7 and 18 are examples 
of this combining technique. 

Amplifying remarks on these programs and subroutines can 
be found in the comments in the computer FORTRAN listings. 

For additional description of plots and computer printed out- 
puts, refer to the appropriate appendices. 
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III. PROGRAM LOGIC AND THEORY 



ECTRACE traces the three-dimensional paths of a bundle of 
rays, one ray at a time. The bundle is divided into fans of 
rays of the same initial elevation angle. A separate ray fan 
is traced for each fixed increment of initial elevation angle 
between the limits ELST and ELEND (input constants) . A posi- 
tive elevation angle describes a ray pointed downward. The 
initial ray headings in each ran fan extend from KDST to HDEND 
(input constants) , measured in degrees clockwise from grid 
north . 

A. BOTTOM AND SOUND FIELD MODELING 

Figure III.l shows the projection of the bottom surface 
onto the horizontal x-y plane. 




FIGURE III.l Bottom Surface Projection 
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The projection of each triangular facet is a right isosce- 
les triangle. The entire surface can be defined by specifying 
the depth values at the vertices of the projection triangles. 
ECTRACE stores these depth values (km) in the two-dimensional 
working array ZB. The parameter DHN (km) represents the spac- 
ing between the vertices in both the X (East) and Y (North) 
directions. Thus, the value of the matrix point ZB(I,J) is 
the depth at the vertex point with X coordinate (I - 1) * 

DHN and Y coordinate (J - 1) * DHN, relative to the ZB origin. 

The matrix values for ZB are read from a source file stored 
on a permanent storage device such as a disk pack. A matrix of 
dimensions 151 by 151 with spacing DHN = km has been used for 
trial runs. 

1 . Bottom Contact Point 

The subroutine CONTAC locates the coordinates and calcu- 
lates the depth of the bottom bounce point. The parameters 
describing locations of the ray segment head and tail are among 
those passed to the main program. In this subroutine, the ray 
segment is treated as unrefracted (a straight line) , a reasonable 
approximation for the small distances involved in determining 
the contact point. 

First the horizontal positions of the ray segment head (x,y) 
called CEE and CNN in CONTAC, and (x x ,y 2 ) and called CRE and 
CRN in CONTAC, are checked to see if they both lie in the same 
projection triangle. If not, an interative routine searches 
for two points on the segment such that the horizontal coordi- 
nates of both points lie inside the same triangular cell and 
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define the portion of the segment that penetrates the bottom 
facet. Then the depth of the ray head z, and the depth of the 
tail z, are computed. 



Figure III. 2. illustrates a profile of a cell bottom with 
a ray intersecting it in the propagation plane. The solution 
of the triad ( X C 'Y C ' Z C ) identifying the point of bottom contact 
proceeds as follows: 

Using subroutine DEEP (described in Section 3) compute and 
H 2 , the bottom depth at {x lf y ) and C X 2 »Y 2 )/ respectively. Then 
compute : 

Az = z - z 
1 2 

where ( X c /Y c ' z c ) are contact point parameters. 

Hz - Hz 

z = 1 2 2 1 

C A _ , TT " 

If AZ ¥■ 0, M 



If AZ = 0, M 

Then x = x + 

c 

y c = y + 

where cf> is the ray heading (angle between the y-z plane and the 
propagation plane) . 

2 . Projection Triangle Vertex Subscripts 

The subroutine IDTSUB calculates the I subscripts, IRT, 
ILT , IVR and the J subscripts JBT , JTP , JVR for the vertices of 



1 2 

- z 

= __C 1_ 

AZ 

= z i - H i 

H, - H 2 

M r sin <p and 

M r cos <}> 
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FIGURE III. 2. Profile of the bottom cell in the propagation 

plane, with a ray contact. 



of the projection triangle which contains the ray's horizontal 
coordinates x and y. The subscript parameters are indicated 
in Figure III. 3. The following suffixes identify the relative 
locations of the vertices of the projection triangle: 

RT - right 



LT - left 
TP - top 
BT - bottom 



■VJZjI'JZ 





1R.T, Iff 



TSLXj J0T 




XK.T, TTP 



IlTjTQT 

FIGURE III. 3. Projection Triangle Vertices 



iv/a, rva 
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3 . Depth Calculation 



Subroutine DEEP calculates the depth in kilometers 
corresponding to the horizontal coordinates x and y. Auto- 
matically, IDTSUB is called to identify the projection tri- 
angle subscripts. The depth z is calculated by solving the 

c 

equation of a plane defined by the depths at the vertex points 
(D j , D 2 , and D^) . Figure III. 4 illustrates the calculations 
made based on the orientation of the projection triangle out- 
lining the bottom facet plane. The quantities x and y repre- 
sent the horizontal position (in units of DHN) of the depth 
calculation point relative to the horizontal position of the 
depth calculation point relative to the horizontal position 
of the depth value D : (the lower left vertex of the projec- 

tion triangle) as calculated in Secion 1. The depth values 
D^ and D^ are normalized by the length of the equilateral 
side 1, in units of DHN. Finally, the depth is solved by 



= x D + 
c x 



y D + D, 
2 c y 1 



B. RAY LOCATION AND DIRECTION PARAMETERS 

A three dimensional left hand coordinate system is used 
(Fig. III. 5.), with the positive z axis pointing in the direc- 
tion of increasing depth and the sea surface level lying in 
the x-y plane. As the ray is traced through three dimensional 
space, the program variables which are used to locate the tail 
of the ray vector are: 

DEPl - depth coordinate . 

CRE - x coordinate in km 
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Or 




FIGURE III. 4. Depth calculation on bottom facet 
CRN - y coordinate in km 

CRE and CRN are components of the horizontal range increment 
DR which forms the propagation plane heading angle PHI with 
the y-axis . 
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Given that the ray vector originates at the point with 
the above coordinates, the direction of the ray vector is 
specified by the following angles: 

TH1 - ray elevation angle in the plane of propagation 
measured from the ray vector's projection on the xy 
plane to the ray vector. A ray pointing toward the 
ocean bottom will have a positive elevation angle. 

Theta takes on values between — tf / 2 and +n/2 (output 
values are in degrees) . 

PHI - ray (propagation plane) heading angle between the 
positive y-axis and the projection of the ray vector on 
the xy plane, measured clockwise from the positive y-axis 
(grid North). PHI takes on values between -t and A. 

Initial values for THl and PHI are calculated from user input 
values (in degrees) , which are converted to radians for compu- 
tational purposes. The unit vectors in the horizontal plane 
are given as UVX = sin (PHI) and UVY = cos (PHI) . 

C. RAY PATH CALCULATIONS 
1 . Path Length 

ECTRACE confines refraction to the vertical plane alone. 
It is recognized that while horizontal sound speed gradients 
exist, their effects are sufficiently small within a well defined 
water mass to be neglected compared to the bathymetric effect. 

However, while the horizontal gradients are likely to 
be slight, important discontinuities in water characteristics 
may occur horizontally across oceanic fronts. For this reason, 
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ECTRACE has been designed to simulate up to three distinct 
water masses and two associated frontal boundaries. 

Within a water mass having uniform characteristics 
along fixed depths a ray path conveniently defines a vertical 
plane within which Snell's Law can be written as 

V = c(z)/cos0(z) (1) 

where V is treated as invariant within the water mass (termed 
the vertex velocity in some texts since it is the speed at the 
depth at which the ray vector becomes horizontal) . The term 
c(z) is the sound speed at depth z and 9 (z) is the elevation 
angle at depth z. 

Equation (1) uses the simplification that the speed of 
propagation and the elevation angle are functions of depth only. 
It is generally possible to define a stratified medium in which 
each layer's gradient is constant. Although such a model re- 
quires a large number of data points to adequately describe a 
water column of complicated sound speed structure, the piece- 
wise linear gradient approximation greatly simplifies the ray 
path descriptions. 

Taking advantage of the approximations, we use a simple 
equation for the sound speed gradient, within a layer, 

g = dc/dz (2) 



Integrating yields 



z = c/g 



(3) 
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where the constant of integration is avoided by arbitrarily 
specifying the origin of the coordinate axes at the depth where 
c = 0. By applying Snell’s Law we obtain 

z =-¥- cos6 ( 4 ) 

and 

dz sin6d0 ( 5 ) 

Specifying r as the horizontal distance axis of the vertical 
plane and applying analytical geometry there results 





dz = tanQdr 


(6) 




dr = — — cos0d6 

g 


(7) 


and 


V . „ 

r = - sinfi 

g 


(8) 


when measuring r from the vertical (z-axis) to the point at 
which 0 = 0 (horizontal projection). 

Squaring (4) and (8) and summing yields 




2 2 , V . 2 

r 2 z 2 = ( g ) 


(9) 


This is the equation of a circle of radius (V/g) 
center is at the origin of the coordinate system 
The equation 


and whose 
just described 




R = - V 

g 


(10) 



describes the radius of curvature of a ray path in a constant 
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sound speed gradient. The minus sign has been chosen to allow 
3 to be positive when the refracted ray is increasing its ele 

vation angls. 

The coordinate system and computer-approximated ray 
path segments are shown in Fig. III. 6. 







FIGURE III. 6. Ray Path Geometry 

It can be shown from Fig. III. 6. that the ray segment 
path length between P and Q is 



A 2, = 0 A0 



UD 



where 



A0 = e 2 - e i 



( 12 ) 



Since the initial 
successive values 



elevation angle of the ray is 
of 0 are found iteratively: 



presumed known, 




(13) 
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npl e geometry also allows for solution of the depth ohange, 
Az = z 2 - Z 1 = P(COS 0 1 - cose 2 ) » (14) 



( 15 ) 



, d likewise the increase in horizontal range, 

Ar = r - r = p(sin9 - sin0 ^ • 

2 1 

«. that in the Fig. ^ ^ 

o the Sign of p must agree with that of .3. Eguation ,10) 

hich gives p the opposite sign of the gradient, satrsfred thr 

eguirement. Only Az is allowed to become negatrve, as w e 

■ay decreases its depth. 

2 . Tra vel Time 

7 rvia tit 7 ) the sound speed relation 

From the diagram (Fig. III./-J 

ship becomes 

df 



c = 3t 
dz 



( 16 ) 

( 17 ) 



dt = 



csin9 
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The ray invariant 



_ 1_ _ cos 9 
3 V c 



( 18 ) 



and the trigonometric identify sin© 
and give dt as a function of depth 



- V 



1 m o 

- cos 9 combine 



dt = 



dz 



( 19 ) 



c(t) \|i - [ac ( z ) ] z 



The time-of-f light of the acoustic ray along the path is found 
by integrating from the initial point to the final point: 



h - h ■ 



dz 



z . 

l 



c(z) \|l - [ac (z) ] 



( 20 ) 



Before integrating it is assumed that the gradient is constant 
along the segment. More generally, ECTRACE assumes a constant 
gradient in the layer bounded by z^ and z ^ . 

Then 



c(z)=c(z)+g(z - z ) 

1 2 1 



( 21 ) 



where z < z £ z 
1 2 



The integral may now be evaluated by introducing a new variable 



W = SlZL such that 



„ _ c(z) 

W = Z - Z ! + — 



( 22 ) 



Then dW = dz and C (W) 



W f 

t- t. = / 

f l w 



= gW. Therefore, 
dw 



i gW \ 



1 - a 2 g 2 w 2 
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( 23 ) 



1 . W [1 + V. 1 - a 2 g 2 Wc ] 

= - In r f_ ' — 1 

g w. [1 + r 1 - a z g z wf "T" 



, W, (1 + sinQ.) 

1 * n [ _J -±- 1 



(24) 



VJ. (1 + sin0 f ) 



Values for N are restricted to being in the same layer and on 
the same side of a turning point in the ray path. 

D. RAY TRACING PROCEDURE 

ECTRACE begins tracing a ray from its given initial coordi- 
nates CRE, CRN, and DEP1. First the water mass is identified 

to determine the local gradient and therefore the ray invarrant 

TiVio rav oath is constructed in the 
and radius of curvature. The ray patn i 

form of arc segments using an increment DTH calculated from a 

chosen arc length DL. Then the new 6 value is calculated and 

Equations (14, and (If, are used to determine the ray segment's 

depth increment and horizontal range increment. The horizontal 

position is updated using the horizontal unit direction vectors 

i i-or) -from PHI the initial heading of the 
(UVX and UVY) calculated from P / 

ray's plane of propagation. 

, , p o D coordinates of the ray head may now 
The interim update 3 

be stated in FORTRAN as 

CEE = CRE + DR * UVX 

(26) 

CNN = CRN + DR * UVY 

(27) 

DEP2 = DEPl and J3Z 

pnts the horizontal range increment Ar and DZ is 
where DR represents the 

the depth incrementAz. 
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E. RAY CONSTRAINT CONDITIONS 



Each time the ray path is incremented, certain constraint 
conditions must be checked. These constraints are bottom re- 
flections and the sound channel limits. 

1 . Bottom Reflections 

After the interim horizontal coordinates CEE and CNN 
are calculated, they are passed to DEEP to find the bottom 
depth WH2 at the ray head. If DEP2 is greater than or equal to 
WH2, a bottom bounce has occurred. In this case the values of 
DEP2, CEE, and CNN are adjusted by linear interpolation and a 
point of bottom contact is established. Subroutine BOANG is 
then called to calculate new values of TH1, PHI, UVY, and the 
ray grazing angle GRAZ. 




FIGURE III. 8. Bottom Bounce Reflection Angles. 

The basic calculations of BOANG involve the vectors: 
v - unit ray vector before bounce 

n — unit vector normal to triangle facet at the bounce point 
v'- unit ray vector after bounce 
The components of v are calculated from TH1 and PHI . The com- 
ponents of n are calculated from the equation of the plane of 
the bottom facet. The new vector v is then calculated from 
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the vector equation, 



v' = v - 2 (v • n) n . (28) 

New values of THl, PHI, UVX and UVY are then calculated from 
the compnents of v' . Finally, the grazing angle GRAZ comes 
from 



sin(GRAZ) = - (v • n) . (29) 

Reflection from a sloping bottom plane facet, as op- 
posed to a surface reflection, forces a reevaluation of the 
Snell's Law constant V, and therefore the radius of curvature 
R. The invariant V remains in effect until the next bottom 
bounce occurs or a new water mass is entered. Thus, after 
summoning BOANG, the ray tracing continues with the new ray 
parameters. The sound channel limits CHS and CHD described 
below are also recalculated at this point since they also 
depend on V, and GRAZ is passed to function FBLOSS to calcu- 
late a bottom loss value. 

2 . Sound Channel Limits 

The Snell's Law equation (1) says that the ray trace 
is constrained to depths where the speed of sound does not 
exceed V. A given ray path may be constrained in the sense 
that it is purely refracted (R) , refracted or reflected (RSR 
or RBR) or purely reflected (BRSR) . The determining factors 
are V and the local sound speed by profile. V is calculated 
from the sound speed at the initial depth and the initial 
elevation angle, and is recalculated after a bottom bounce or 

upon entering a new water mass. 
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CHS and CHD are the refractive sound channel limits 
at the shallow and deep extents respectively. Since the sound 
speed profile is approximated by piecewise continuous linear 
gradients, the water column is treated as a column of indexed 
layers, each characterized by upper and lower (shallow and 
deep) depth boundaries and a constant sound speed gradient. 

The limiting depths CHS and CHD are found by scanning the 
water mass layers both above and below the current ray loca- 
tion until a boundary sound speed value is found which exceeds 
V. The limiting depth is then calculated from 

z, . = z + (V - c)/g (30) 

lim 

where z and c are evaluated at the top of the layer and g is 
the gradient of the limiting layer. 

For cases where g = 0 or where V is not exceeded by 
any sound speed in the local profile, CHS and CHD are set to 
the physical limits of the water column. 

F. RAY PROPAGATION TERMINATION CONDITIONS 

Tests are made during each ray segment iteration for events 
which terminate the ray trace. These events are: 

• Escape from grid - when the ray departs the finite x and 
y constraints of the working depth array (ZB) . 

• Depth boundaries - when the depth detected at the ray is 
less than ten wavelengths or greater than 10 km. 

• Bottom loss - when the attenuation due to bottom reflec- 
tion exceeds a specified input amount in dB. 
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• Error at bottom contact - when defective program code 
or anomalous depth data prevent the interpolative rou- 
tine CONTAC from recovering the bottom depth within a 
reasonable number (50) of iterative attempts. 

• Trapping - when it becomes apparent that the ray will 
not contact the bottom, rendering further three-dimen- 
sional tracing unnecessary. 

• Bottom bounce gate - when the total number of bottom 
reflections exceeds a user specified amount. 

The cause of ray termination is specified in the printed 
ray history, giving the ray parameters at the time of 
termination . 

G. IMPROVEMENTS MADE OVER THE ORIGINAL NRL VERSION 

The ray tracing algorithm originally existed as a program 
called ABOUNC at the Naval Research Laboratory. ECTRACE dif- 
fers significantly from ABOUNC in the following areas: 

1 . Ray Path Development 

ABOUNC calculates ray parameters through right tri- 
angle solutions and Snell's Law, but develops the ray path 
along fixed depth. Since this prevents rays from becoming 
horizontal, a reversal is forced when the magnitude of the 
elevation angle goes below a gate value. In contrast, ECTRACE 
develops the ray path along arc segments and radii of curva- 
ture, fixing only the maximum length of the segments while 
interrupting the increment as necessary when a segment reaches 
a reversal point or layer boundary, or to keep the increment 
within a maximum amount . 
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2 . Computer Word Storage 



ABOUNC floating point variables are nearly all FORTRAN 
double precision. The largest single burden on computer core 
is the depth matrix, read from a magnetic tape in card image 
format, each card record containing eight items (fields) ten 
bytes long. This format permitted a depth range of ±9999 km 
with resolution of ±1 cm. By reducing the depth range to ±999 
km and resolution to ±10 cm, two unnecessary bytes are removed 
from each field, enabling a fit of ten items per record and a 
permanent storage reduction of twenty percent. In addition, 
ECTRACE reads the matrix into a single precision array to re- 
duce the core requirements for the data by fifty percent. 

3 . Grid File Manipulation 

ABOUNC reads the contents of the entire depth matrix 
data file into array whose dimensions must be exactly equal to 
those of the matrix. This requirement is not inconvenient on 
a large computer with FORTRAN dialect which permits objects 
time dimensioning. ECTRACE is written in FORTRAN G for an IBM 
Series/360 Model 67 machine whose limitations required some 
fundamental modifications of the program and its job control. 

To maintain flexibility and reduce run time, ECTRACE 
was reduced to a package of subprograms which have been pre- 
compiled and stored on a permanent device in load module form. 
The user must supply a calling program which serves primarily 
to establish the dimensions of the depth array. The user's 
selection of array dimensions are made on the basis of core 
economy and area of interest rather than the dimensions of the 
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input matrix. ECTRACE, through its internal subroutine GRDSCT, 
extracts a rectangular working subregion from any part of the 
input matrix. Trial runs of ECTRACE used a working array of 
dimensions 151 by 151 representing a sector of 150 km on a side 
with a core requirement of less than 250K, while the input ma- 
trix was 369 by 443. Without GRDSCT, ECTRACE would have re- 
quired an allocation of 800K bytes using the same input file. 

4 . Treatment of Small Elevation Angles 

As stated previously, ABOUNC forced a reversal of the 
ray as the elevation angle of the leading segment approached 
a gate value close to zero. Since the ECTRACE tracing scheme 
permits a full range of elevation angles, early trial runs 
revealed an inherent error phenomenon as 6 approached zero in 
a steep gradient. Specifically, a A 0 on the order of 3 de- 
grees caused by a small radius of curvature or along arc seg- 
ment traced a ray to its vertex depth in error as much as 50 
meters from the Snell's Law prediction. Since the depth change 
was calculated from the difference in elevation angle cosines, 
the error was found to arise from the machine inability to 
retain difference precision for any pair of numbers close to 
integer values (round off error) . Since the difference rather 
than the actual cosine values is needed for this calculation, 

0 2 a 4 

the small angle series approximation cos 1 - - . * * 

was used to formulate the equation 

cos 6 - cos 6 = %(0 2 -0 2 ) - — (0 4 -0**). ( 31) 

1 2 2 1 2 “* 2 1 

When restricted to small angles, this approximation yields 
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greater accuracy than the straightforward machine calculation 
of the left hand side, since the precision retained by the 
registers is very high. 

5 . Travel Time 

ABOUNC calculates the travel time of an acoustic sig- 
nal along a ray segment by dividing the linearly-approximated 
path length by the mean sound speed. ECTRACE accounts for 
path curvature and depth dependence of the sound speed through 
integral evaluation, discussed in detail in section C.2. 

6 . Track Events 

Layer transitions, reflections, and turnarounds are 
events which may require a change in ray parameters. ABOUNC 
tests for these events in algorithmic sequence, makes the 
necessary changes in response to the first event detected 
and continues the trace from the current position of the ray 
vector head. ECTRACE retraces a ray segment to the event 
depth and tests for all other events before beginning the next 
segment from the transition point. This logic eliminates the 
possibility of a long ray segment crossing more than one event 
depth with only one event detected. 
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IV. USER INSTRUCTIONS (ECTRACE) 



A single step ECTRACE job traces a bundle of rays from a 
single source in the form of ray fans. Each ray fan is a bun- 
dle of rays of common initial elevation angle. All ray fans 
are bounded between to initial heading values measured from 
grid North (grid North equals true North at the grid center) . 

During the following discussion the reader may wish to 
refer to Fortran listing in Appendix M. The user instructions 
for the other programs (GENBOT, G3DP , SYNGEN, GRDCHK , and 
ECCOM) are explained in their individual appendices (User In- 
structions and Output Description) . 

A. JOB CONTROL CARDS (JCL) 

The following JCL cards and parameters are required for 
an ECTRACE run; 

• JOB card - The CPU time parameter should allow 60 sec- 
ccnds for every seven rays to be traced. 

• EXEC card - This card must specify the FORTCLGW proce- 
dure (Appendix 0) . 250K bytes of core is required for 

the standard 22801 (151 by 151) element depth matrix 
portion plus 4K for every additional 1000 elements. 

The computer systems job class definitions should be 
considered before deciding upon the size of the working 
matrix and the number of rays. For example, 

EXEC FORTCLGW, REGION. GO=25 OK 
is a class C job at NPS when limited to five minutes of 
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time, which is adequate for tracing 40 rays. 

• Calling program - Following the 

/ /FORT . SYS IN DD * 

card, the calling program must contain cards 
COMMON/ DIMS /KMX , KMY and 
CALL TRACER (ZB) . 

ZB is a REAL * 4 array dimensioned exactly by KMX, KMY, and 
these values (KMX, KMY) are initialized to the values 
of the DIMENSION statement using assignment statements, 
as in the following example. 

DIMENSION ZB ( 151 , 151) 

COMMON /DIMS/ KMX, KMY 

# 

« 

KMX=151 

KMY=151 

CALL TRACER (ZB) 

STOP 

END 

KMX and KMY are INTEGER * 4 and should be less than or 
equal to the dimensions of the data matrix from which 
the depth values will be taken. 

• LINK cards - These cards point to the ECTRACE load 
module. They must include all members of the module 
for which the user does not supply a substitute. An 
example which includes all members is: 

//LINK . USDD UNIT=333 0 , VOL=SER=DISKOl , DSN=127 0 . ECTRACE , 

// DISP=SHR, LABEL= ( , , , IN) 

/ /LINK . SYSIN DD * 
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INCLUDE USDD (TRACER, GRDSCT , IDPROF , NUPROF , CHNLIM) 

INCLUDE USDD (BGNPLT , SSPPLT , BDTPLT , RNGPLT , T2DPLT , ENDPLT ) 
INCLUDE USDD ( IDTSUB , DEEP , CONTAC , BOANG , ANGPRT , FBLOSS ) 
ENTRY MAIN 
/* 

♦ GO cards - ECTRACE uses data set reference number one 
for input of the depth matrix. For example, to point 

to a depth matrix residing in a data set called S9999.GRID 
on DISK03, the following DD statement is used: 

//GO . FT01F001 DD UNIT=330, VOL=SER=DISK03 ,DSN=S9999 . 

GRID , DISP=SHR, LABEL= ( , , ,IN) 

If a depth matrix from an outside source is used, it 
must be in the format specified in the Output Grid Data 
Format section of Appendix A. If a punched card output 
is not desired to make an ECCOM run, include the card 
//GO.FT07F001 DD DUMMY 

ahead of the GO.FTOlFOOl DD statement shown above. 

B . DATA CARDS 

Following the 

//GO.SYSIN DD * 

card, the user selected data must be provided in the following 
format. FORTRAN field descriptors are in parentheses. The 
FORTRAN variable names are listed for reference. 

* Card 1 ( 20A4 ) 

OTIT - title of the run 

• Card 2 (958.2) 

XO - source position for x coordinate (kilometers) . 
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YO 



ZO 



XREC 

YREC 

ZREC 



FREQ 



FDL 



source position y coordinate (kilometers) . 
source depth (meters) . A source depth which 
exceeds the water depth will be decreased to one 
meter above the bottom. Source coordinates in a 
nonpropagating location will also be adjusted, 
receiver position x coordinate (kilometers) . 
receiver position y coordinate (kilometers) . 
receiver depth (meters) . The receiver may be 
positioned at the source if calcuations are not 
desired. Excessive depth is automatically 
corrected . 

frequency (Hz) . A ray is terminated when the 
water depth becomes less than ten wavelengths 
during propagation. 

(optional) segment length (meters) . The parameter 
affects the smoothness of the vertical plot. 
Smaller values increase smoothness and run time. 
Large values create the risk of missing a bottom 
peak. The default value is one wavelength. 



♦ Card 3 (6F8.2) 



XMIN - initial x value of horizontal plot (kilometers) . 
YMIN - initial y value of horizontal plot (kilometers) . 
The above values are used for drawing the axis of the 
horizontal plot and for determining the section of the 
input matrix to be read into a working depth matrix (ZB) . 
The values for KMX and KMY, recorded cell resolution 
(DHN) , and the limits of the input matrix ultimately 
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determine the axis ranges. The user should be familiar 
with these parameters when choosing the origin of the 
plot. 

ZST - initial depth value for vertical plot (kilometers) 
ZEND - final depth value for vertical plot (kilometers) 
ARST - initial range value for vertical plot (kilometers) 
AREND - final range value for vertical plot (kilometers) 

The above values are used to draw the axes of the vertical 
plot. The initial values may be negative. The program 
adjusts the final values if necessary for plot esthetics. 

• Card 4 (6F8.2) 

HOST - heading of the first ray in each ray fan (degrees) 
HDEND - heading of the last ray in each ray fan (degrees) 
ELST - elevation angle of the first ray fan (degrees) 
ELEND - elevation angle of the last ray fan (degrees) 
Headings are measured from grid north. Elevations are 
measured from the horizontal plane, positive downward 
(toward increasing depth) . 

HORLIM - maximum ducting or channeling range (kilometers) . 

This is a gate value, like BLMAX CARD 2, which 
terminates the trace of a ray when exceeded. 
HORLIM serves to screen out those rays which 
do not appear to contact the bottom. 

SMNGL - (optional) small angle value for alternative 
cosine difference calculation (degrees) . A 
value of less than three degrees is recommended. 
See the discussion on "small angle approximations" 
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in Chapter III. The default value of two 
degrees is set when left blank. 

Card 5 (615) 



NHD 


number of rays per ray fan (tl) 


NEL 


number of ray fans per bundle (>1) 


NSTOP - 


maximum number of bottom bounces per ray (>1) 



The program restricts NHD to 30 rays, with a maximum 
total number of 100 rays produced, each with a maximum 
number of 50 bottom bounces or 200 reversals. 



INCPLT 


- number of ray segments defining one plot seg- 
ment (>1) . With sufficiently small ray seg- 
ments it is not necessary to use each iteration 
for plot definition to achieve a smooth plot, 
so this parameter is provided to enable the 
user to coordinate FDL and INCPLT for run 
time economy. 


NLEN 


- length of the longest horizontal plot axis 
(<19 inches) . This value controls the physi- 
cal size of the plot. A value less than four 
will supress the plot (resulting in an ABEND 
code in the HASP log) , but it will not affect 
the calculations or printed and punched out- 
put. A value greater than 15 requires addi- 
tional JCL for strip plotting, since the 
total plot width is NLEN plus six inches. 
Refer to Ref. 7 for strip plotting. 
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NSSP - number of sound speed profiles (water masses 
provided as data ( 1 < NSSP < 3) . 

• Card 6 (A4,I3) 

OMAS (1) - name of first water mass 

NVAL (1) - number of (depth, sound speed) pairs defining 
the sound speed profile (1 £ NVAL £ 50) . 

• Card 7 - (through card 16 if necessary) (10F8.2) 

These are the data pairs (depth, sound speed) , in meters 
an m/sec., from the sea surface to the last data point 
in the profile. The gradients will be calculated linearly 
between these points, with the final gradient below the 
last point set to 0.01668 m/sec. by the program. 

• Next card (4F8.2) 

If the previous data is to be followed by the data for 
another water mass, this card must be included to define 
the frontal boundary between the two. The order in 
which the water masses appear as inptut determines their 
positions relative to the fronts. Cards which describe 
the next water mass, follow the frontal definition card 
and are formatted as described cards 6 and beyond. The 
last data card is the last string of (depth, sound speed) 
pairs of the last water mass; refer to Appendix M (sec- 
ond page) for an example of three water masses. The 
following variables define the boundary between each 
water mass: 

FENTX (1) - x coordinate of front beginning point (km) 

A 

FENTY (1) - y coordinate of front beginning point (km) 
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FENTX (1) - x coordinate of front end point (km) 

FEXTY (1) - y coordinate of front end point (km) 

The examples shown in Fig. IV. 1 identify the water 
mass boundaries for various input axis values in the 
form (FENTX, FENTY) , (FENTX, FEXTY) for each water mass 
boundary. The first example a. represents the data 
from an original ECTRACE run, refer to the second page 
of Appendix M. Exambles b. and d. show the priority 
of the water masses. The first water mass defined, 
Northwestern Mohns Ridge (NWMR) , dominates the second. 
Southeastern Mohns Ridge (SEMR) , and the third, isove- 
locity (ISOV) , where they overlap. 
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a. 




NWMR 

(0,117) (88,150) 
SEMR 

(38,0) (150,50) 
ISOV 




NWMR 

(0,108) (150,108) 
SEMR 

(0,46) (150,46) 
ISOV 



b. 




NWMR 

(21,0) (108,150) 
SEMR 



(0,118) 

ISOV 


(124,0) 




d. 


NWMR 


f ^ 1 11 

/ 
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// 

1 


\ 
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\ 

\ 


1 


/ 

SEMR / 

i 


ISOV 


I 

i 

... 1 





31 73 9fl 150 



NWMR 

(0,83) (98,150) 
SEMR 

(31,0) (73,150) 
ISOV 



FIGURE IV. 1. Examples of water mass boundaries for various values of 
(FENTX,FENTY) and (FEXTX,FEXTY) . 
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V. INTERPRETATION OF ECTRACE OUTPUT 



The following discussion refers to an ECTRACE run titled 
"Trial 67***Synthetic Basin***" (Figure 5) , which is a trace 
conducted in a synthetic parabolic basin (Figure 6) . A por- 
tion of the printed output is shown in Tables II - V. In this 
ECTRACE run, eight ray fans were selected to be projected be- 
tween -30 and +30 degrees elevation. Each ray fan was com- 
posed of only one ray, thus the initial heading for each ray 
was 010 degrees. 

A. PRINTED OUTPUT 

Refer to Table II through V (Appendix G) as examples of 
the following discussion. 

1 . Parameter Tableau (Table II ) 

After processing and adjusting the input parameters 
just described, ECTRACE prints a final tableau of these parame 
ters plus some additional parameters calculated from the origi 
nal input such as layer gradients, the total number of rays 
to be traced, and the spacing between rays. In addition, the 
results of the depth matrix-grid load is summarized, revealing 
data from the header record of the grid file. 

2 . Ray History Tableau (Table III) 

During the trace of each individual ray, this tableau 
prints a line of ray parameters at each end point and reversal 
The following is a description of one line of parameters. 

* TYPE identifies the nature of the reversal, and is 
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abbreviated INIT for initial position, BOTT for bottom 
bounce, SURF for surface reflection, OVER for turnover, 
and STOP for terminated position. 

• X, Y, and DEPTH are the coordinates in kilometers of the 
reversal and end point. 

• GRID HEADING is the ray heading measured clockwise from 
grid North at the completion of the reversal. 

• ELEV is the ray's elevation angle measured downward 
(positive toward increasing depth) from the horizontal 
zero degree reference) . 

• BPUG TO REC is the bearing from the reversal point to 
the receiver position. 

• DIST TO REC is the slant range distance from the reversal 
point to the receiver. 

• WTR MASS is the water mass identifier indicating the 
sound speed profile applied in the present calculations . 

• LYR is the constant gradient layer number corresponding 
to the order input as profile data (i.e., layer 1 is the 
surface layer and indicates that the surface gradient 
applies . 

• DEPTH LIMITS are the physical and refractive sound chan- 
nel limits of the ray, determined by its invariant. If 
MIN is 0.0, an upward traveling ray will reflect off 
the surface before reaching its upper vertex (turnover) 
depth. If the MAX depth is greater than or equal to 10km, 
a line of asterisks is printed since the ray could never 
attain its turnunder depth before reaching the bottom. 
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Otherwise, these limits are the ray's vertex depths, 
and if the reversal type is OVER or UNDR, the appropri- 
ate limit should exactly match the depth value on the 
same line. 

• C is the sound speed at the reversal point. This is the 
ray invariant used for tracing. 

• GRAZ ANGLE is the angle in degrees between the ray vector 
and the line of intersection of the bottom facet plane 
with the vertical plane. This is the angle used to 
calcualte bottom loss and the new ray direction vector. 

• BTM ANGLE is the angle between the bottom facet plane 
normal and the vertical (z-axis) , in degrees. 

• BTM LOSS is the accumulated intensity loss in dB after 
the bottom reflection. When this value equals or exceeds 
BLMAX, the ray's trace will terminate. 

• LPS is the number of steps required by the CONTAC sub- 
routine to locate the point of bottom contact. A long 
ray segment length may require more steps than a shorter 
one. A maximum number of 50 steps is permitted before 
the bottom facet is assumed anomalous and the ray's 
trace is terminated. 

• Prior to the final line of an individual ray's trace 
history, is a summary of plot definition parameters plus 
data on the ray's CPA to the receiver, and reason for 
ray trace termination. The CPA calcualtions are accurate 
to within one ray segment length. 
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3. Ray Fan Summary Tableau (Table IV) 



At the completion of the trace of one complete ray 
fan, a table is printed listing each ray's elevation angle 
and depth in meters before each bottom bounce. All zero values 
at the bottom of each ray’s column indicate blank data. 

4. CPA Summary Tableau (Table V) 

When all ray traces have been completed, a table is 
printed listing each ray's heading, time, and distance at CPA 
to the receiver. 

B. PLOTTED OUTPUT 

ECTRACE produces a single montage of up to five separate 
plots (Fig,. il3 , the total dimensions of which depend upon the 
user's selection for the maximum size in inches of the longest 
axis of the horizontal plot (NLEN) . The following table lists 
some sample values for NLEN and the resulting dimensions of 
the complete plot. 



NLEN 


Maximum 
(in. ) Height 


Maximum 

Width 


Remarks 


5 


7 


8.33 


minimum 


10 


14 


16.67 




15 


21 


25 


maximum size without strip 
plotting [5] 


19 


26.6 


31.67 


maximum allowable 



Components of the ECTRACE plot follow. 

1. Sound Speed Profiles (Figure 13) 

A plot of sound speed versus depth is made for each 
water mass defined. 



55 



2. Vertical Plot (Figs. 15 and 16) 



Most two-dimensional ray trace programs produce plots 
of ray paths on graphs of depth versus range. The vertical 
plot is similar except that all ray paths are projected from 
three dimensions onto the vertical plane. The plane is oriented 
along the mean ray fan heading and includes its profile of the 
sea bed. In ECTRACE trial 65B (Figure 13) the ray fans each 
contain two rays whose initial headings vary from 360 to 020 
degrees as indicated by 10.000 +/- 10.000 on the plot. The 
bathymetry profile is made along the heading of 010 degrees. 

On a plot, the heading value (HEADINGS) always indicates the 
direction of the bathymetry slice. The (+/-) value signifies 
the initial heading of the first and last rays and that addi- 
tional rays are evenly spaced in the interval. 

Since the rays are actually traced in three dimensions, 
bottom reflections occurring outside the vertical plot plane 
may appear to be taking place above or below the bottom sur- 
face instead of at the boundary (Fig. 15). Only those bottom 
contact points which occur in the projection plane must lie 
on the seabed profile. 

The symbols listed on the plot (Fig. 5) are used to 
identify bottom contact locations in the vertical and horizon- 
tal plot, and are unique for each ray fan except that up to 
thirteen symbols are used before repetition, as shown in 
Table VI. 

3. Horizontal Plot (Fig. 17) 

This plot is the most important plot product of 
ECTRACE and is the plot of the ray paths on a horizontal plane 
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(topview) . This plot displays the X-Y positions of each 
bottom contact point using the ray's (ray fan) identification 
symbol and connects them with straight lines, revealing the 
horizontal deflection effect. Each ray's plot begins at the 
source, ends up at the last bottom contact point in the grid 
and is labeled with the ray number (Fig. 17) . The ray num- 
ber can be used to locate its printed and punched output 
history. 



The axes are labeled with respect to the origin of 
the input matrix. The geographic center position and coor- 
dinates of the input matrix, if included in the bounds of 
the working matrix (FORTRAN array ZB) , is also plotted for 
reference. It is important to note that this plot, without 
the bottom contact symbols, may combine the results of other 
ECTRACE runs by combining the optional punched card outputs 
as input to ECCOM (Figs. 7 and 18). 



C. PUNCHED OUTPUT 

The punched card output for a single ECTRACE job consists 
of the following cards: 

CARD DATA CONTENTS 

1 The runtitle 

2-3 The grid title 

4 Grid reference coordinates and scale 

factor 

5 Initial elevation, initial heading, 

number of points (BB) and ray num- 
ber of the first ray 

6+ X-Y coordinate pairs of plot points 

of the first ray 

(Data for additional rays is repeated as 5 and 6+ above) 

Last 9999., to flag the end of data 
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These cards may be stacked as a member with those of 
several ECTRACE jobs as input to a single ECCOM job. For 
ECCOM to be meaningful, cards 2 through 4 of each stack mem- 
ber should be identical, and the user should verify this 
himself . 

Punched card output may be quite voluminous and should 
not be produced unless required. The following card, inserted 
before the first GO card (JCL) , will suppress the punched 
output : 

//GO.FTQ7F001 DD DUMMY 
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VI . 



SUMMARY 



A. POTENTIAL APPLICATIONS 

ECTRACE, in its present form, has many potential 
applications : 

• The investigation of asympotic deflection angles in 
ocean areas where bathymetry data are well-documented, 
allows optimum hydrophone positioning. For example, 
Ref. 6 cites cases where seemingly attractive receiver 
positions would be affected by shadowing in the hori- 
zontal plane. 

• Adaptation of ECTRACE ray tracing methods could be 
incorporated into existing programs which presently 
rely upon an assumption of radial symmetry of ray 
propagation about a source or receiver (no horizontal 
deflection) . 

• The horizontal deflection effect may lend additional 
insight into the study of travel times and intensities 
of acoustic signals, enabling more accurate source 
fixing information. 

B. IMPROVEMENT 

There are many areas of improvement which would increase 
the utility of ECTRACE and its supporting computer programs 
for the above applications: 

• A better routine to convert bathymetry contour data 
into a matrix of depth values is essential for modeling 
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regions for which only contour data exists. Most known 
interpolation algorithms perform best on data which are 
randomly or regularly spaced rather than in the form of 
contours. For example RAIN2 [5] would probably perform 
better with the data used to make the contours than it 
does with the contours themselves. While the contour 
data are useful for cartographic and geologic purposes, 
bathymetry data in a less processed form are significantly 
better sources for constructing computer models of the 
sea bed. 

ECTRACE in its present form makes no calculations of 
intensities or their losses due to spreading, attenua- 
tion, or scattering, thus preventing its use as a 
propagation loss model. 

An improved bottom loss function is recommended over 
the supplied FORTRAN function sabprogram FBLOSS. The 
user may perform this substitution easily by deleting 
FBLOSS from the INCLUDE step in the LINK statement of 
the ECTRACE job control and supplying another function 
subprogram (added at the end of the calling program) 
that accepts a double-precision argument for grazing 
angle and is also called FBLOSS. 

The assumption of specular reflection neglects the 
phenomenon of transmission into the sediment and its 
effect on the location of the reflected ray due to sub- 
bottom boundaries. It would be interesting to incorporate 
in an improved version of ECTRACE. 
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APPENDIX A 

GENBOT USER INSTRUCTIONS AND OUTPUT DESCRIPTION 

GENBOT is used in the construction of a depth grid file from bathymetry 
contour files. In this discussion, the term "grid" refers to a matrix of 
depth values whose row and column dimensions represent data point spacing 
in terms of a desired resolution parameter DHN. 

1. Selection of Bathymetry Grid Regions 

GENBOT will create a bathymetry grid file for any region of the geoid 
for which contour data is supplied in the format stated in section A-4. 

The user must select the latitude and longitude of the center of the 
region of interest, a radius in kilometers and the cell resolution in 
kilometers. GENBOT will attempt to create a square matrix whose semi- 
axes represent the selected radius and whose elements represent depth 
values in kilometers. 

A 30-km data margin, outside the calculated axis limits, is required 
to begin interpolation. Should this requirement not be met, GENBOT will 
symmetrically reduce the length of one or both axes until the resulting 
rectangular matrix plus margins with the desired center position intact 
is totally encompassed within the data region. The center coordinates, 
area radius, and cell resolution must be selected so as to avoid generating 
a grid which is impracticably small. In particular, a generated grid may 
not overlap a pole, and one with a pole on a margin boundary must be sup- 
plied data from eighteen contour files to cover 180 degrees of longitude. 

Contour data currently available in usable format cover 40 degrees 
of longitude of the North Atlantic Ocean, specifically the Norwegian and 
Greenland Seas, from longitudes 20W to 20E and latitudes 60N to 90N 
(Fig. 8). The Naval Research Laboratory, Acoustics Division, Arctic Section, 
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has additional files which cover the entire North Atlantic land and sea 
regions which may be utilized by GENBOT after modification. However, 
GENBOT need not be the exclusive means of grid instruction. The only 
requirement is that the grid exists as a file in the format specified. 

The user may wish to draw on other souces of bathymetry data, such as 
SYNBAPS [8], and tailor other available resources to construct a usable 
grid. 

2 . Grid Construction Techniques 

Once supplied with the user’s selection for the grid title, center 
coordinates, cell resolution and area radius (TTL, CNLAT , CNLON, DHN, 
ARAD), GENBOT begins by defining a north-oriented square region whose 
boundaries are at a distance ARAD from the center. The subscripts of 
the two dimension matrix depicting this region will represent equal units 
of DHN spacing from the origin along the x and y axes of the square using 
a cartesian coordinate system. Thus, while all distances calculated 
during grid generations are geodetic distances from the center, point 
positions are made relative to the grid origin or southwest corner. This 
technique results in a plane projection which is distortion free along 
straight tracks which pass through the center. Subroutine CORNER then 
calculates and prints a table of the initial cartesian and geographic 
coordinates of the grid corners. The cartesian coordinates will always 
range from (0,0) to their maxium values, initially double the value of 
ARAD. The geographic coordinates should reflect some meridian conver- 
gency. This table will be repeated for the final grid, showing the 
adjustments made, if any, to meet the data requirements. 

Many geographic coordinate pairs in an input file may be well outside 
the area relevant to the grid, so GENBOT calculates gate values called 
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SLAT, ELAT, SLON, and ELON for data discrimination. SLAT and ELAT are 
always one degree beyond the minimum and maximum latitude values calcu- 
lated by CORNER, but SLON and ELON may extend beyond the minimum and 
maximum longitudes by 1, 10, or 20 degrees, depending upon the maximum 
latitude. 

The grid region and its data environs recognized by GENBOT at this 
stage is shown in Fig. A.l. All points within these environs are sub- 
ject for consideration during interpolation of depth values in the grid 




FIGURE A.l. Geographic region used for grid construction. 
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Data files used as input for GENBOT are digitized represenations of 
an NRL bathymetry chart [4] and were recorded on Mobil Oil Company Tape 
NAV-78 for use in a computer program to draw portions of the chart in the 
form of bathymetric contour lines and their depth labels as they exist 
on the chart. The records are in the card image format and consist of 
three types: 

• NAME - depth value coordinates 

• CON - contour line header label 

• Chart coordinates of points along the contour line 

In the original form the numerical data on the tape consisted of lati- 
tude/longitude pairs, providing no means of determining the depth repre- 
sented by the contour lines without visual inspection of the computer plot. 
The tape used in this research is similar in format to NAV-78 except that 
the contour line header records include depth data (provided by NRL) , and 
only the regions shown in Fig. A. 2 and listed in Table A. I, are covered. 

The numerals of the data set names were chosen to reflect the original 
file number on the source tape, NAV-78. 

Using the magnetic tape files selected by the user, GENBOT begins 
reading the digitized contour data, stopping after each record to dis- 
criminate the points, calculate geodetic distances, transform the data 
to x, y, z coordinates, and load the coordinates into three axis vectors, 

AX, AY, and AZ (which will be used later for interpolation) while simul- 
taneously keeping a running plot of the contour lines. 

After all input files have been processed, the vectors AX and AY are 
sorted into ascending order by the SSP3 subroutine RAIN1 [5]. The mini- 
mum and maximum values are tested to determine if the 30-km outside data 
margin is included in the range. If not, axis reduction begins until 
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FIGURE A. 2. North Atlantic regions, with their magnetic tape file numbers 
and data set names, available for use by GENBOT. All are 10° 
by 10°. Other regions available on the same tape are listed 
in Table A. I . 
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TABLE A. I 

FILE ONE OF CONTOUR DATA RECORDED ON A MAGNETIC TAPE 

DATA: BATHYMETRIC CONTOURS OF THE NORWEGIAN AND GREENLAND SEAS FROM 

LAT 6 ON TO 90N and LONG 20W to 20E. 

SOURCE: ADAPTED FROM MOBILE OIL CO. TAPE NO. NAV-78 AND NAVAL RESEARCH 

LABORATORY DATA. 

STRUCTURE: EACH FILE COVERS A 10 X 10 DEGREE REGION. THE FOLLOWING LIST 

IDENTIFIES THESE FILES. THE FIRST VALUE IS THE FILE NUMBER, FOLLOWED BY 
THE DSN, THEN THE LAT/LONG PAIR OF THE SOUTHWEST CORNER (LABOR, LONOR) OF 
THE REGION, AND FINALLY THE NUMBER OF CARD IMAGES CONTAINED IN THE FILE. 
THE TAPE IS STANDARD LABELED, DENSITY 1600 BPI , LRECL=80,BLKSIZE=800. 



2 


C0NF13 


+80-20 


256 


5 


C0NF26 


+80-20 


4343 


8 


CONF34 


+70+00 


12443 


11 


CONF35 


+70-10 


9428 


3 


CONF24 


+60+10 


569 


6 


CONF27 


+60-10 


6001 


9 


CONF37 


+80+00 


491 


12 


CONF36 


+80+10 


409 


4 


CONF25 


+60+00 


3026 


7 


CONF30 


+70+10 


857 


10 


CONF06 


+70-20 


4454 


13 


CONF38 


+80-10 


972 
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either the margin requirement is met or goes below a 7.5 km semi-axis 
length, the latter case preventing grid construction although a contour 
plot may still be produced. Subroutine CORNER is then called again to 
print the corner coordinate table. If no adjustments were necessary, the 
second table will be identical to the first. 

The interpolation process begins at the southwest corner of the grid 
and proceeds in a columnar direction (i.e., completing the range of x^ 
values for each y value) via successive calls to the SSP3 subroutine 
RAIN2. The accuracy of the interpolated values is proportional to the 
density of data in the given area of interest. Every ten values inter- 
polated are written as a record on the output file. The process continues 
until the total number of matrix elements has been reached, at which time 
the output file is closed and the program is terminated. 

3. User Instructions 

a. JCL and Program Deck. 

The program section of a sample deck is printed below. The TIME 
and REGION. GO parameters were chosen in anticipation of a maximum number 
of 5000 data points in the area of interest, a square region of 210 km 
on a side (150 km plus margin) , The program listings included in this 
appendix are part of this example. 

//(Standard job card) ,TIME=20 

// EXEC F0RTCLGW, REGION. GO=250K (See Appendix 0.) 

// FORT. SYS IN DD * 

(GENBOT, CORNER, GE0PLT , GE0DST listings) 

(RAIN1 , RAIN 3 , RAIN 2 listings) 

/* 

The following DD statements identify the magnetic tape files 
representing the ten by ten degree regions into which the subroutine 
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extends. In the example, the area chosen is centered at coordinates 
70°-04 'N/15°-51 ' (70.07059/15.85302) with a radius of 75 km. This area 
requires two data files as input, CONF24 (file 3) and CONF30 (file 7). 

Since tape data is sequential, the order of the DD cards should be in the 
order of the file numbers (appearing in the LABEL parameter). 

//GO . FT01F001 DD UNIT=3400-3, VOL=SER=NPS,DISP= (OLD, PASS)-, 

// LABEL=(3,SL, ,IN) ,DSN=C0NF30 

//GO . FT01F002 UNIT=AFF=FT01F001 , V0L= ( , RETAIN , REF=* . GO . FT01F001) , 

// DISP= (OLD, PASS) ,LABEL= (7 , SL , ,IN) ,DSN=C0NF30 
More information on magnetic tape usage at NPS may be found in Ref. 10. 

Next must follow the DD statements specifying the data set and 
device for permanent storage of the generated grid. The permanent device 
may be another magnetic tape, a disk, or punched cards. In all cases, 
the size of the output data set is important, and the DCB parameters 
should be considered in advance if the grid size is different from the 
example. GENB0T produces card image out put, meaning the logical record 
length is 80 bytes. 

• Punched cards . If the grid is to be stored on punched cards, the 
user need only change the data set reference number on the appro- 
priate WRITE statements in GENBOT from a two to a seven. The num- 
ber of cards punched will equal a 3+[(row dimension) (column 
dimension) /10] rounded up to an integer value. This represents 
three header cards followed by ten depth values per card in columnar 
order, as required by ECTRACE. Thus if the input grid radius was 
75 km and input cell resolution 1 km, the maximum parameter calculated 
by GENBOT will be: 

XMAX YMAX IMAX JMAX 

150 150 151 151 
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resulting in a matrix of 22801 elements and a punched output of 
2284 cards. 



• Disk . Grid storage on a disk device is the most convenient form for 
multiple ECTRACE runs over the same sea bed. The following DD state- 
ment stores the example grid on a 3330 disk drive device at NPS as a 
data set named S999.GRID. 

//GO . FT02F001 DD UNIT=3330,V0L=SER=DISK03,DISP=(,KEEP) , 

// DCB= (FECFM=FB , LRECL=80 , BLKSIZE=800) ,SPACE= (TRK, 20) , 

// DSN=S9999 . GRID 

Blocked at 800 bytes, the output data requires 14 tracks per 10 
records (card images). More information on the DCB and DSN parameters 
may be found in Ref. 11. When the generated grid is no longer needed, 
and if the user wishes to generate a new grid, or if space has already 
been allocated on the device, the DISP parameter should be SHR in 
place of (KEEP) and the DCB and SPACE specifications may be omitted. 
Disk data sets are periodically purged by center operations personnel, 
but obsolete data set space should be cleared or reallocated by the 
user . 

b. Input Data. 

The following input data are user selected values and follow 
the //GO.SYSIN DD* card. 

Card Variable FORTRAN Field Descriptors Remarks 

1-2 TTL 6A8/6A8 Title 

3 CNLAT, 4F10.5 center latitude, 



CNLON , 
DHN, ARAD 



center longitude 
cell resolution, 
radius 



4 



NFILES, LEN 



215 



Number of input 
files, plot 
length 
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The magnitude of LEN determines the plot size (£19 inches) . If 
LEN is negative, no grid will be constructed. If its magnitude is less 
than four, no plot will be produced. 

c. Printed Output 

GENBOT printed output supplies the following: 

. A repeat of the user input data. 

• A table of the grid corner coordinates and the depth matrix 
dimensions using the input parameters. 

• A summary of results for each contour line, including its depth, 
and the title from the header label, the number of points on the 
line (NPTS) and the number of points used from the line (KPTS) . 

• A corner coordinate table of the final grid. 

• A summary of the interpolation results. A line is printed for 
each interpolated point whose accuracy may be questioned. The 
meaning of the IER parameter may be found in the RAIN1 listing 
in Appendix I. 

d. Plotted Output 

The contour lines of the input file are plotted with the local 
parallels and meridians overlaid. The resultant subregion represented by 
the output grid is outlined by a rectangle with labeled axes. The axis 
values are referenced to the desired grid origin and not the resultant 
grid origin, the latter being at (0,0) always. 

e. Programming Notes on Contour Data Interpolation 

Subroutine RAIN 2 , the interpolation routine used by GENBOT, is 
aptly suited to interpolate points that are more or less randomly spread 
in both horizontal directions. The nature of contour data, however, is 
typically a collection of closely spaced depth values along widely spaced 
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contour (compare Figures 111,1. and A.I.). This type of data arrangement 
frequently results in a colinearity problem for RAIN 2 , which can be resolved 
only by searching ever widening data cells for a non-colinear point, result- 
ing in a degradation of accuracy. In addition, the routine is highly sensi- 
tive to contour line kinks which result in right-angle deviations in the 
interpolated contour lines (which may be plotted by GRDCHK) . 

Correcting the deficiencies of an interpolating routine is beyond 
the scope of this research, but the Naval Research Laboratory, Arctic Sec- 
tion, is currently investigating methods for generating grid-matrices from 
contour data. Users are encouraged to modify the code used in GENBOT as 
desired to produce results which more closely reflect the actual data 
supplied . 

4. Input Contour Data Format 

As stated previously, some contour data already exists in a format 
usable by GENBOT, but if data is supplied by another source, it must be 
arranged in a file in the manner and format specified below. All records 
are 80 byte card images. 

The first record of the file contains data on the file itself. GENBOT 
uses a 415 field to represent the origin latitude/longitude pair, the num- 
ber of contour lines in the file, and the number of card images it may 
expect to read. 

Contour line data exists as separate sequences of cards and may 
appear in any part of the file as long as each line or line segment is 
in an unbroken sequence with a header card. The header card must contain 
the line’s depth in hundreds of meters, right justified in a 14 field in 
columns 41-44. The cards following the header card contain the latitude/ 
longitude pairs in degrees and decimals of degrees of the points along 
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the contour line in the form of four pairs of coordinates (in FORTRAN 
fields 8F10.5) per card. South and west values are negative. A value 
greater than 90.0 in magnitude in any latitude field is interpreted as the 
end of data for the contour line, and GENBOT expects the next record to 
be either another contour line header label, the end of the file, or a 
depth value position record. 

Depth value position records are optional and are used during the 
chart-drawing phase of the program to place alphanumeric symbols on the 
chart product. They may exist any place after the first record in the 
file except imbedded in a contour line sequence. Columns 4-8 must con- 
tain "NAME" in an A4 field, columns 19-38 contain the latitude and longi- 
tude of the symbols in a 2F10.5 field and the symbols themselves must be 
in columns 46-47. These records are normally used to draw depth values 
in hundreds of meters on the contour lines. 

If the user wishes to use GENBOT on externally supplied contour data, 
it it recommended that a dump be performed on one of the smaller files 
(such as C0NF13) of the GENBOT data tape to verify the required format 
for adaptation. 

5 . Output Grid Data Format 

GENBOT constructs the output file in accordance with GRDSCT require- 
ments. Since GRDSCT is the subroutine used by TRACER (of ECTRACE), 

G3DP and GRDCHK, any externally supplied grid data must all be of the 
same structure as outlined below in order to be used by those programs. 
Again, all records are card images. 

a. The first two records contain the input grid title in the 1QA8 
fields. 

b. The third record contains the input latitude and lontitude of the 
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grid center, the point spacing, and the number of rows and columns of the 
matrix representing the grid. The record format is (3F10 . 5 , 215) . 

c. The fourth and remaining records contain ten depth values per 
card, in kilometers, as 10F8.4 fields, beginning with column one (y = 0) 
and proceeding all x-values before going to the next column, in accordance 
with the FORTRAN convention of columnar storage of matrices. 
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APPENDIX B 

SYNGEN DESCRIPTION AND USER INSTRUCTIONS 



1. Program Description 

a. SYNGEN is a collection of simple functions which may be selectively 
used to generate and record a matrix of depth values representing a desired 
synthetic bottom configuration. Six bottom cases are provided (see below), 
but the program may be easily modified by the user to produce a case not 
described here. A program listing is provided in Appendix H. 

The user selected parameters are: 

• CNLAT - any desired center latitude value. 

• CNLON - any desired center longitude value. These parameters 

are requested in compliance with' the required output file 
format . 

• DMIN - the minimum depth value (D^) to be attained within a 
75 km radius from the grid center. 

• DMAX - the maximum depth value (D^) within the same radius. 

• ITYPE - an integer code for the bottom type desired. 

b. Table B.l lists the ITYPE code, the bottom type, the functions 

used to generate the depth matrix, and the matrix dimensions. All grids 

represent a region covering a 150 km in both the x and y directions, but 

the matrix size will be 11 by 11 or 151 by 151, depending upon the bottom 

type (DHN values are 15 km and 1 km, respectively). In all cases, 

AD = D - D . 
l o 

c. The depth file provided by SYNGEN conforms to the format required 
by ECTRACE, GRDCHK, and G3DP . When one of these programs uses a SYNGEN 
product, the array dimension in the program T s source code must be checked 
to conform to the dimensions of the generated matrix. 
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TABLE B.I 

SYNTHETIC SEA BED FUNCTIONS 



ITYPE Type Function 



Parameters 



12 Wedge 



13 Ridge 



H = D - Y 

l 



H = D +|y6(x) | 
o 

H = D elsewhere 
o 



Y = AD/135 

0 £ x £ 135 

Y = AD/15 
6(x) = x - 75 

60 £ x £ 90 



14 Trough 



H = D -|y<$(x) | 6 = AD/15 

H = D i elsewhwere 6(x) = x - 75 

60 < x < 90 



15 Conical H = D + yr y = AD/75 

Seamount 0 t ; 

H = D elsewhere r =\I(x - 75) 2 + (y -75) 

o 1 

16 Sinusoidal H = D y [l+cos{6(x) /R} ] y = Ad/ 2 
Undulation 

S(x) = 2(x - 75) 

R = 15 



17 Parabolic H = D 
basin 



yr 



Y = AD/ (75) 2 
r = nJ(x - 75) 2 + (y - 






Matrix 

Dimensions 

11 X 11 
11 X 11 

11 X 11 

151 X 151 
151 X 151 

151 X 151 
75) 2 
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2. User Instructions 



THE JCL of a sample deck are arranged as follows: 

//(Standard JOB card),TIME=2 
// EXEC FORTCLG , REGION .GO=250K 
/ /FORT . SYSIN DD * 

(SYNGEN source code) 

/* 

/ /GO.FT02F001 DD UNIT=3330,VOL=SER=DISK03 , 

// DSN=S9999 . SYNBED , DCB= (RECFM=FB ,LRECL=80 ,BLKSIZE=800) , 

// DSP= (NEW, KEEP ) , SPACE= (TRK, 20) 

//GO. SYSIN DD * 

(user selected data) 

/* 

Notice that the output data set is directed to DISK03. The DSN is the 
user's selection for the data set name in accordance with NPS system re- 
quirements. A 151 x 151 depth matrix requires twenty tracks of disk 
space on a 3330 device. 

Only one data card is needed to supply the program with the user's 
parameter selections. Its field descriptor is (4F10.5,I5) and contains 
the values reflected in the same order as listed in paragraph l.a. 



76 



APPENDIX C 

G3DP DESCRIPTION AND USER INSTRUCTIONS 
1 . Program Description 

a. G3DP uses the NPS library subroutine PLT3D1 [Ref. 12] to draw a 
perspective surface plot of a subset of a depth matrix. This plot is a 
normal projection of the surface using a focal length and line-of-sight 
calculated from the user's selections for the viewpoint position angles 
DEGUP and DEGCW. As the variable names imply, the viewer is elevated 
DEGUP degrees from a reference depth and rotated DEGCW clockwise from a 
south-to-north viewing aspect. For example, the standard working matrix 
(dimensioned 151 by 151) is best viewed with DEGUP and DEGCW values of 

30 degrees and 45 degrees respectively, providing the same viewing aspect 
demonstrated in all the the G3DP products printed in this report. Using 
a DEGUP value of zero results in an isometric rather than a perspective 
projection. The DEGCW selection should be between zero and 89 degrees 
(south-north to west-east viewing aspects) . 

b. As in ECTRACE and GRDCHK, G3DP requires source code modification 
by the user to tailor the program to the desired dimensions of the work- 
ing matrix. Instructions in this regard are included as comments in the 
FORTRAN listing (Appendix K) . 

c. An important characteristic of G3DP is that it uses only the 
odd-numbered rows and columns of the working matrix to produce the sur- 
face plot. This reduction in resolution is required to reduce core re- 
quirements and unclutter the graphic output caused by typically large 
(greater than 75 by 75) working matrices. 
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2. User Instructions 



The JCL of a sample program deck are arranged as follows: 

//(Standard JOB card) ,TIME=2 

//EXEC FORTCLGW, REGION. GO=300K (Appendix 0.) 

//FORT . SYSIN DD* 

(G3DP source code) 

/* 

/ /LINK. USDD UNIT-3330, V0L=SER=DISK01, DSN- 1570 . ECTRACE, 

// DISP=SHR,LABEL=( , , , IN) 

//LINK. SYSIN DD* 

INCLUDE USDD (GRDSCT) 

ENTRY MAIN 

/* 

//G0.FT01F 

//GO. SYSIN DD* 

(user-selected data) 

/* 

The LINK step provides access to GRDSCT for loading the working 
matrix from the input file. The input data set is the same example used 
in ECTRACE and GRDCHK. The card input data are user selected and con- 
tained on one card, as follows (FORTRAN field descriptor 5F8.3) 

• XST - working grid matrix x-origin relative to that of the 
input grid (km) . 

• YST - working grid-matrix y-origin relative to that of the input 
grid (km) . 

• DEGUP - viewpoint elevation angle (0 <, DEGUP < 89 degrees). 
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• DEGCW - viewpoint rotation angle from south to north 
(0 £ DEGCW < 89 degrees) . 

• LW - width of the plot (5 < LW < 18 inches). 



79 



APPENDIX D 

GRDCHK-USER INSTRUCTIONS AND OUTPUT DESCRIPTION 



1 . Purpose and Output Description 

Cases may arise where the accuracy of depth values in a generated ma- 
trix may be in question because of data sparsity, interpolation deficiencies 
or other reasons. Program GRDCHK assists in evaluating the generated grid 
by printing and plotting selected columns of the matrix for comparison with 
that of the source contours. The contour plot may also be used on a light 
table with an ECTRACE or ECCOM horizontal-ray trace plot (Figs. 5, 6, and 
7). 

The NPS library routine CONTUR is called for producing the contour 
plot. Subroutine GEOPLR is a modified version of GEOPLT to accommodate 
the rotated plot produced by CONTUR. . 

2 . User Instructions 

To aid in this discussion refer to a computer listing of GRDCHK code 
in Appendix J. A listing of the CONTUR subroutine may be found in Ref. 6. 

The following example applies to a standard 151 by 151 working grid. 
Below is a listing of the JCL cards and the position of the program deck. 

//(Standard job card), TIME=5 

// EXEC FORTCGLW, REGION. GO=350K 

//FORT. SYS IN DD * 

(GRDCHK, CORNER, GEOPLT, GEODST source decks) 

/* 

The following DD statements link to a member of ECTRACE for allocating 
a working subregion of the input. 

//LINK. USDD DD UNIT=3330,VOL=SER=DISK01, 
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/ / DSN= 15 7 0. ECTRACE , DISP=SHR , LABEL= ( , , , IN) 

/ /LINK. SYSIN DD * 

INCLUDE USDD (GRDSCT) 

ENTRY MAIN 

/* 

Next follows the DD statement identifying the depth matrix input file. 

The same example is used for the GENBOT output file. 

/ /GO.FTOFOOl DD DSN=S9999 .GRID,UNIT=3330,VOL=SER=DISK03, 

// DISP=SHR, LABEL- (,,, IN) 

Finally the user input parameters are included. The card 
//GO. SYSIN DD * 

is followed by: 

Card Parameter Names FORTRAN field description 

1 NL, LEN 215 

2 ZST, ZEND, XMIN, 4F6.1 

YMIN 

3 IPS, IPE, JPS, JPE 415 

where 

• NL is the number of contour levels desired. The shallowest and 
deepest levels, if known, should be read into ZST and ZEND. If 
the boundary levels are not known, NL should be made negative, 
and subroutine CONTUR will calculate the boundary and contour 
levels. If NL is greater in magnitude than 100, no contour plot 
will be produced. 

* LEN is the length of the contour plot's longest side in inches. 

If LEN is negative, a cartesian grid will be drawn on the con- 
tour plot instead of geographic coordinate lines. If LEN is 
less than 4, no contour plot will be produced. The maximum 
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magnitude of LEN is 19 inches. 

• ZST AND ZEND are estimates of the smallest and largest values in 
the working matrix (ZB) . 

• XMIN and YMIN establish the origin of the working matrix relative 
to that of the input matrix. The program assignments for KMX and 
KMY, the recorded cell spacing DHN, and the dimensions of the work- 
ing matrix determine the subregion size parameters in the same man- 
ner as they do in ECTRACE, since subroutine GRDSCT is linked to 
this program. In effect, the user must tailor certain portions 

of the program cards to compensate for the FORTRAN G inability 
to perform object time dimensioning. 

• IPS, IPE, JPS and JPE identify the beginning and end rows and 
columns of the working depth matrix to be singled out for inspec- 
tion. Thus, if ZB is dimensioned (150,150) and the above values 
are 1, 150, 74, and 78 respectively, GRDCHK will print the depth 
values of these 750 elements of the ZB matrix and plot five 
columnar profiles (Fig. 4). If this record is omitted, only 

the contour plot will be produced. 
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